####Model results####
##a) minority-majority units (table X15, models A45-A48)
min_main_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant_g == 1))
min_main_m1g_t1_sc_cse <- data.frame(cluster.se(min_main_m1g_t1_sc, as.factor(subset(main_group,  subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
min_main_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant_g == 0))
min_main_m1g_t1_sf_cse <- data.frame(cluster.se(min_main_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
min_main_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant == 1))
min_main_m1d_t1_sb_cse <- data.frame(cluster.se(min_main_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant == 1)$cowcode))[, 2])
min_main_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant == 1))
min_main_m2d_t1_sb_cse <- data.frame(cluster.se(min_main_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(min_main_m1g_t1_sc, min_main_m1g_t1_sf, min_main_m1d_t1_sb, min_main_m2d_t1_sb, se=c(min_main_m1g_t1_sc_cse, min_main_m1g_t1_sf_cse, min_main_m1d_t1_sb_cse, min_main_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A45)", "Min. (model A46)","Maj./min. dyad (model A47)","Maj./min. dyad (model A48)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X15. Additional results: minority-majority units only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex15.txt")
stargazer(min_main_m1g_t1_sc, min_main_m1g_t1_sf, min_main_m1d_t1_sb, min_main_m2d_t1_sb, se=c(min_main_m1g_t1_sc_cse, min_main_m1g_t1_sf_cse, min_main_m1d_t1_sb_cse, min_main_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A45)", "Min. (model A46)","Maj./min. dyad (model A47)","Maj./min. dyad (model A48)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X15. Additional results: minority-majority units only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex15.html")
##b) majority-majority units (table X16, models A49-A51)
maj_main_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant_g == 0))
maj_main_m1g_t1_sf_cse <- data.frame(cluster.se(maj_main_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
maj_main_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant == 1))
maj_main_m1d_t1_sb_cse <- data.frame(cluster.se(maj_main_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
maj_main_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant == 1))
maj_main_m2d_t1_sb_cse <- data.frame(cluster.se(maj_main_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(maj_main_m1g_t1_sf, maj_main_m1d_t1_sb, maj_main_m2d_t1_sb, se=c(maj_main_m1g_t1_sf_cse, maj_main_m1d_t1_sb_cse, maj_main_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Min. (model A49)","Maj./min. dyad (model A50)","Maj./min. dyad (model A51)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X16. Additional results: majority-majority units only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex16.txt")
stargazer(maj_main_m1g_t1_sf, maj_main_m1d_t1_sb, maj_main_m2d_t1_sb, se=c(maj_main_m1g_t1_sf_cse, maj_main_m1d_t1_sb_cse, maj_main_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Min. (model A49)","Maj./min. dyad (model A50)","Maj./min. dyad (model A51)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X16. Additional results: majority-majority units only.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex16.html")
##c) all violence incidents (table X17, models A52-A55)
allinc_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, int_grp_rel_dominant_g == 1))
allinc_m1g_t1_sc_cse <- data.frame(cluster.se(allinc_m1g_t1_sc, as.factor(subset(main_group, int_grp_rel_dominant_g == 1)$cowcode))[,2])
allinc_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, int_grp_rel_dominant_g == 0))
allinc_m1g_t1_sf_cse <- data.frame(cluster.se(allinc_m1g_t1_sf, as.factor(subset(main_group, int_grp_rel_dominant_g == 0)$cowcode))[,2])
allinc_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, int_grp_rel_dominant == 1))
allinc_m1d_t1_sb_cse <- data.frame(cluster.se(allinc_m1d_t1_sb, as.factor(subset(main_dyad, int_grp_rel_dominant == 1)$cowcode))[, 2])
allinc_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, int_grp_rel_dominant == 1))
allinc_m2d_t1_sb_cse <- data.frame(cluster.se(allinc_m2d_t1_sb, as.factor(subset(main_dyad, int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(allinc_m1g_t1_sc, allinc_m1g_t1_sf, allinc_m1d_t1_sb, allinc_m2d_t1_sb, se=c(allinc_m1g_t1_sc_cse, allinc_m1g_t1_sf_cse, allinc_m1d_t1_sb_cse, allinc_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A52)", "Min. (model A53)","Maj./min. dyad (model A54)","Maj./min. dyad (model A55)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X17. Additional results: including group/dyad-years even if the group/dyad has been involved in civil/communal violence in the previous year.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex17.txt")
stargazer(allinc_m1g_t1_sc, allinc_m1g_t1_sf, allinc_m1d_t1_sb, allinc_m2d_t1_sb, se=c(allinc_m1g_t1_sc_cse, allinc_m1g_t1_sf_cse, allinc_m1d_t1_sb_cse, allinc_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A52)", "Min. (model A53)","Maj./min. dyad (model A54)","Maj./min. dyad (model A55)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X17. Additional results: including group/dyad-years even if the group/dyad has been involved in civil/communal violence in the previous year.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex17.html")
##d) only groups in core area [population share in unit > 0, i.e. int_grp_rel > 0] (table X18, models A56-A59)
only0l_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, int_grp_rel > 0 & subset_analysis == 1 & int_grp_rel_dominant_g == 1))
only0l_m1g_t1_sc_cse <- data.frame(cluster.se(only0l_m1g_t1_sc, as.factor(subset(main_group, int_grp_rel > 0 & subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
only0l_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, int_grp_rel > 0 & subset_analysis == 1 & int_grp_rel_dominant_g == 0))
only0l_m1g_t1_sf_cse <- data.frame(cluster.se(only0l_m1g_t1_sf, as.factor(subset(main_group, int_grp_rel > 0 & subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
only0l_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, int_grp_rel1 > 0 & int_grp_rel2 > 0 & subset_analysis == 1 & int_grp_rel_dominant == 1))
only0l_m1d_t1_sb_cse <- data.frame(cluster.se(only0l_m1d_t1_sb, as.factor(subset(main_dyad, int_grp_rel1 > 0 & int_grp_rel2 > 0 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
only0l_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, int_grp_rel1 > 0 & int_grp_rel2 > 0 & subset_analysis == 1 & int_grp_rel_dominant == 1))
only0l_m2d_t1_sb_cse <- data.frame(cluster.se(only0l_m2d_t1_sb, as.factor(subset(main_dyad, int_grp_rel1 > 0 & int_grp_rel2 > 0 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(only0l_m1g_t1_sc, only0l_m1g_t1_sf, only0l_m1d_t1_sb, only0l_m2d_t1_sb, se=c(only0l_m1g_t1_sc_cse, only0l_m1g_t1_sf_cse, only0l_m1d_t1_sb_cse, only0l_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A56)", "Min. (model A57)","Maj./min. dyad (model A58)","Maj./min. dyad (model A59)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X18. Additional results: only groups in 'core' settlement area.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex18.txt")
stargazer(only0l_m1g_t1_sc, only0l_m1g_t1_sf, only0l_m1d_t1_sb, only0l_m2d_t1_sb, se=c(only0l_m1g_t1_sc_cse, only0l_m1g_t1_sf_cse, only0l_m1d_t1_sb_cse, only0l_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A56)", "Min. (model A57)","Maj./min. dyad (model A58)","Maj./min. dyad (model A59)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X18. Additional results: only groups in 'core' settlement area.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex18.html")
##e) only groups that are considered politically relevant throughout by EPR (table X19, models A60-A63)
always_relevant_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, relevant_throughout == 1 & subset_analysis == 1 & int_grp_rel_dominant_g == 1))
always_relevant_m1g_t1_sc_cse <- data.frame(cluster.se(always_relevant_m1g_t1_sc, as.factor(subset(main_group, relevant_throughout == 1 & subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
always_relevant_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m1g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, relevant_throughout == 1 & subset_analysis == 1 & int_grp_rel_dominant_g == 0))
always_relevant_m1g_t1_sf_cse <- data.frame(cluster.se(always_relevant_m1g_t1_sf, as.factor(subset(main_group, relevant_throughout == 1 & subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
always_relevant_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, relevant_throughout1 == 1 & relevant_throughout2 == 1 & subset_analysis == 1 & int_grp_rel_dominant == 1))
always_relevant_m1d_t1_sb_cse <- data.frame(cluster.se(always_relevant_m1d_t1_sb, as.factor(subset(main_dyad, relevant_throughout1 == 1 & relevant_throughout2 == 1 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
always_relevant_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, relevant_throughout1 == 1 & relevant_throughout2 == 1 & subset_analysis == 1 & int_grp_rel_dominant == 1))
always_relevant_m2d_t1_sb_cse <- data.frame(cluster.se(always_relevant_m2d_t1_sb, as.factor(subset(main_dyad, relevant_throughout1 == 1 & relevant_throughout2 == 1 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(always_relevant_m1g_t1_sc, always_relevant_m1g_t1_sf, always_relevant_m1d_t1_sb, always_relevant_m2d_t1_sb, se=c(always_relevant_m1g_t1_sc_cse, always_relevant_m1g_t1_sf_cse, always_relevant_m1d_t1_sb_cse, always_relevant_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A60)", "Min. (model A61)","Maj./min. dyad (model A62)","Maj./min. dyad (model A63)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X19. Additional results: only groups that EPR codes as 'politically relevant' throughout each country's existence.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex19.txt")
stargazer(always_relevant_m1g_t1_sc, always_relevant_m1g_t1_sf, always_relevant_m1d_t1_sb, always_relevant_m2d_t1_sb, se=c(always_relevant_m1g_t1_sc_cse, always_relevant_m1g_t1_sf_cse, always_relevant_m1d_t1_sb_cse, always_relevant_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A60)", "Min. (model A61)","Maj./min. dyad (model A62)","Maj./min. dyad (model A63)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X19. Additional results: only groups that EPR codes as 'politically relevant' throughout each country's existence.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex19.html")

####Figure A13####
##a) minority-majority units (civil)
me_min_main_m1g_t1_sc <- margins_summary(min_main_m1g_t1_sc, variables = "sa_territory_t", vcov = cluster.vcov(min_main_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant_g == 1)$cowcode)))
me_min_main_m1g_t1_sc$term <- "second-order maj."
me_min_main_m1g_t1_sf <- margins_summary(min_main_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(min_main_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant_g == 0)$cowcode)))
me_min_main_m1g_t1_sf$term <- "second-order min."
me_main_all <- rbind.fill(me_min_main_m1g_t1_sc, me_min_main_m1g_t1_sf)
me_main_all$term <- as.factor(me_main_all$term)
me_main_all$term = factor(me_main_all$term,levels(me_main_all$term)[c(2,1)])
me_main_all$model <- "minority-majority units"
me_main_all$estimate <- me_main_all$AME
me_main_all$conf.low <- me_main_all$lower
me_main_all$conf.high <- me_main_all$upper
##b) minority-majority units (communal)
me_min_main_m1d_t1_sb <- margins_summary(min_main_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(min_main_m1d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant == 1)$cowcode)))
me_min_main_m1d_t1_sb$term <- "all"
me_min_main_m1d_t1_sb$model <- "minority-majority units"
me_min_main_m2d_t1_sb <- margins_summary(min_main_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(min_main_m2d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 0 & int_grp_rel_dominant == 1)$cowcode)))
me_min_main_m2d_t1_sb$term <- ifelse(me_min_main_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_min_main_m2d_t1_sb$model <- "minority-majority units"
me_main_alld <- rbind.fill(me_min_main_m1d_t1_sb, me_min_main_m2d_t1_sb)
me_main_alld <- subset(me_main_alld, term != "other")
me_main_alld$term <- as.factor(me_main_alld$term)
me_main_alld$estimate <- me_main_alld$AME
me_main_alld$conf.low <- me_main_alld$lower
me_main_alld$conf.high <- me_main_alld$upper
##c) majority-majority units units (civil)
me_maj_main_m1g_t1_sf <- margins_summary(maj_main_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(maj_main_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_maj_main_m1g_t1_sf$term <- "second-order min."
me_maj_main_m1g_t1_sf$term <- as.factor(me_maj_main_m1g_t1_sf$term)
me_maj_main_m1g_t1_sf$model <- "majority-majority units"
me_maj_main_m1g_t1_sf$estimate <- me_maj_main_m1g_t1_sf$AME
me_maj_main_m1g_t1_sf$conf.low <- me_maj_main_m1g_t1_sf$lower
me_maj_main_m1g_t1_sf$conf.high <- me_maj_main_m1g_t1_sf$upper
##d) majority-majority units units (communal)
me_maj_main_m1d_t1_sb <- margins_summary(maj_main_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(maj_main_m1d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_maj_main_m1d_t1_sb$term <- "all"
me_maj_main_m1d_t1_sb$model <- "majority-majority units"
me_maj_main_m2d_t1_sb <- margins_summary(maj_main_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(maj_main_m2d_t1_sb, as.integer(subset(main_dyad, subset_analysis == 1 & igrd_state_control == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_maj_main_m2d_t1_sb$term <- ifelse(me_maj_main_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_maj_main_m2d_t1_sb$model <- "majority-majority units"
me_maj_main_alld <- rbind.fill(me_maj_main_m1d_t1_sb, me_maj_main_m2d_t1_sb)
me_maj_main_alld <- subset(me_maj_main_alld, term != "other")
me_maj_main_alld$term <- as.factor(me_maj_main_alld$term)
me_maj_main_alld$estimate <- me_maj_main_alld$AME
me_maj_main_alld$conf.low <- me_maj_main_alld$lower
me_maj_main_alld$conf.high <- me_maj_main_alld$upper
##e) all violence incidents (civil)
me_allinc_m1g_t1_sc <- margins_summary(allinc_m1g_t1_sc, variables = "sa_territory_t", vcov = cluster.vcov(allinc_m1g_t1_sc, as.integer(subset(main_group, int_grp_rel_dominant_g == 1)$cowcode)))
me_allinc_m1g_t1_sc$term <- "second-order maj."
me_allinc_m1g_t1_sf <- margins_summary(allinc_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(allinc_m1g_t1_sf, as.integer(subset(main_group, int_grp_rel_dominant_g == 0)$cowcode)))
me_allinc_m1g_t1_sf$term <- "second-order min."
me_allinc <- rbind.fill(me_allinc_m1g_t1_sc, me_allinc_m1g_t1_sf)
me_allinc$term <- as.factor(me_allinc$term)
me_allinc$term = factor(me_allinc$term,levels(me_allinc$term)[c(2,1)])
me_allinc$model <- "all violence incidents"
me_allinc$estimate <- me_allinc$AME
me_allinc$conf.low <- me_allinc$lower
me_allinc$conf.high <- me_allinc$upper
##f) all violence incidents (communal)
me_allinc_m1d_t1_sb <- margins_summary(allinc_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(allinc_m1d_t1_sb, as.integer(subset(main_dyad, int_grp_rel_dominant == 1)$cowcode)))
me_allinc_m1d_t1_sb$term <- "all"
me_allinc_m1d_t1_sb$model <- "all violence incidents"
me_allinc_m2d_t1_sb <- margins_summary(allinc_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(allinc_m2d_t1_sb, as.integer(subset(main_dyad, int_grp_rel_dominant == 1)$cowcode)))
me_allinc_m2d_t1_sb$term <- ifelse(me_allinc_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_allinc_m2d_t1_sb$model <- "all violence incidents"
me_allincd <- rbind.fill(me_allinc_m1d_t1_sb, me_allinc_m2d_t1_sb)
me_allincd <- subset(me_allincd, term != "other")
me_allincd$term <- as.factor(me_allincd$term)
me_allincd$estimate <- me_allincd$AME
me_allincd$conf.low <- me_allincd$lower
me_allincd$conf.high <- me_allincd$upper
##g) only groups in core area (civil)
me_only0l_m1g_t1_sc <- margins_summary(only0l_m1g_t1_sc, variables = "sa_territory_t", vcov = cluster.vcov(only0l_m1g_t1_sc, as.integer(subset(main_group, int_grp_rel > 0 & subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_only0l_m1g_t1_sc$term <- "second-order maj."
me_only0l_m1g_t1_sf <- margins_summary(only0l_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(only0l_m1g_t1_sf, as.integer(subset(main_group, int_grp_rel > 0 & subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_only0l_m1g_t1_sf$term <- "second-order min."
me_only0l <- rbind.fill(me_only0l_m1g_t1_sc, me_only0l_m1g_t1_sf)
me_only0l$term <- as.factor(me_only0l$term)
me_only0l$term = factor(me_only0l$term,levels(me_only0l$term)[c(2,1)])
me_only0l$model <- "only groups in 'core' area"
me_only0l$estimate <- me_only0l$AME
me_only0l$conf.low <- me_only0l$lower
me_only0l$conf.high <- me_only0l$upper
##h) only groups in core area (communal)
me_only0l_m1d_t1_sb <- margins_summary(only0l_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(only0l_m1d_t1_sb, as.integer(subset(main_dyad, int_grp_rel1 > 0 & int_grp_rel2 > 0 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_only0l_m1d_t1_sb$term <- "all"
me_only0l_m1d_t1_sb$model <- "only groups in 'core' area"
me_only0l_m2d_t1_sb <- margins_summary(only0l_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(only0l_m2d_t1_sb, as.integer(subset(main_dyad, int_grp_rel1 > 0 & int_grp_rel2 > 0 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_only0l_m2d_t1_sb$term <- ifelse(me_only0l_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_only0l_m2d_t1_sb$model <- "only groups in 'core' area"
me_only0ld <- rbind.fill(me_only0l_m1d_t1_sb, me_only0l_m2d_t1_sb)
me_only0ld <- subset(me_only0ld, term != "other")
me_only0ld$term <- as.factor(me_only0ld$term)
me_only0ld$estimate <- me_only0ld$AME
me_only0ld$conf.low <- me_only0ld$lower
me_only0ld$conf.high <- me_only0ld$upper
##i) only groups that are always politically relevant (civil)
me_always_relevant_m1g_t1_sc <- margins_summary(always_relevant_m1g_t1_sc, variables = "sa_territory_t", vcov = cluster.vcov(always_relevant_m1g_t1_sc, as.integer(subset(main_group, relevant_throughout == 1 & subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_always_relevant_m1g_t1_sc$term <- "second-order maj."
me_always_relevant_m1g_t1_sf <- margins_summary(always_relevant_m1g_t1_sf, variables = "sa_territory_t", vcov = cluster.vcov(always_relevant_m1g_t1_sf, as.integer(subset(main_group, relevant_throughout == 1 & subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_always_relevant_m1g_t1_sf$term <- "second-order min."
me_always_relevant <- rbind.fill(me_always_relevant_m1g_t1_sc, me_always_relevant_m1g_t1_sf)
me_always_relevant$term <- as.factor(me_always_relevant$term)
me_always_relevant$term = factor(me_always_relevant$term,levels(me_always_relevant$term)[c(2,1)])
me_always_relevant$model <- "only groups that are\nalways 'relevant'"
me_always_relevant$estimate <- me_always_relevant$AME
me_always_relevant$conf.low <- me_always_relevant$lower
me_always_relevant$conf.high <- me_always_relevant$upper
##j) only groups that are always politically relevant (communal)
me_always_relevant_m1d_t1_sb <- margins_summary(always_relevant_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(always_relevant_m1d_t1_sb, as.integer(subset(main_dyad, relevant_throughout1 == 1 & relevant_throughout2 == 1 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_always_relevant_m1d_t1_sb$term <- "all"
me_always_relevant_m1d_t1_sb$model <- "only groups that are\nalways 'relevant'"
me_always_relevant_m2d_t1_sb <- margins_summary(always_relevant_m2d_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(always_relevant_m2d_t1_sb, as.integer(subset(main_dyad, relevant_throughout1 == 1 & relevant_throughout2 == 1 & subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_always_relevant_m2d_t1_sb$term <- ifelse(me_always_relevant_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_always_relevant_m2d_t1_sb$model <- "only groups that are\nalways 'relevant'"
me_always_relevantd <- rbind.fill(me_always_relevant_m1d_t1_sb, me_always_relevant_m2d_t1_sb)
me_always_relevantd <- subset(me_always_relevantd, term != "other")
me_always_relevantd$term <- as.factor(me_always_relevantd$term)
me_always_relevantd$estimate <- me_always_relevantd$AME
me_always_relevantd$conf.low <- me_always_relevantd$lower
me_always_relevantd$conf.high <- me_always_relevantd$upper
##k) combined (civil)
figurea13_plot_data_g <- rbind.fill(me_main_all,me_maj_main_m1g_t1_sf,me_allinc,me_only0l,me_always_relevant)
figurea13_plot_data_g$model <- as.factor(figurea13_plot_data_g$model)
figurea13_plot_data_g$term <- as.factor(figurea13_plot_data_g$term)
figurea13_plot_data_g$term = factor(figurea13_plot_data_g$term,levels(figurea13_plot_data_g$term)[c(2,1)])
figurea13_plot_g <- dwplot(figurea13_plot_data_g,
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(colour = model,shape =model)), 
                      whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("a) civil violence") +
  scale_color_manual(name = "Coefficient for:", values = c("#7fc97f","#beaed4","#fdc086","#386cb0","#f0027f"))+ 
  scale_linetype_manual(name = "Coefficient for:", values = c(5,4,3,2,1))+ 
  scale_shape_manual(name = "Coefficient for:", values = c(15,16,17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of civil violence") + ylab("group type")+ #+
  guides(shape = guide_legend(nrow=2,"model",reverse=TRUE), colour = guide_legend(nrow=2,"model",reverse=TRUE),linetype = guide_legend(nrow=2,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format())
##l) combined (communal)
figurea13_plot_data_d <- rbind.fill(me_main_alld,me_maj_main_alld,me_allincd,me_only0ld,me_always_relevantd)
figurea13_plot_data_d$model <- as.factor(figurea13_plot_data_d$model)
figurea13_plot_d <- dwplot(figurea13_plot_data_d,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("b) communal violence maj./min. dyad") +
  scale_color_manual(name = "Coefficient for:", values = c("#7fc97f","#beaed4","#fdc086","#386cb0","#f0027f"))+ scale_linetype_manual(name = "Coefficient for:", values = c(5,4,3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(15,16,17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type")+ #+
  guides(shape = guide_legend(nrow=2,"model",reverse=TRUE), colour = guide_legend(nrow=2,"model",reverse=TRUE),linetype = guide_legend(nrow=2,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format(), breaks = c(0, 0.005, 0.01))
##m) combined
figurea13 <- grid_arrange_shared_legend(figurea13_plot_g, figurea13_plot_d, ncol=2, nrow=1)
ggsave(figurea13, file='../figures/figurea13.pdf', width = 14, height = 7, units="cm",dpi=1000)